Code
knitr::opts_chunk$set(
message = FALSE
)Arbuscular mycorrhizal fungal communities in French vineyard: a signature of terroir whatever the farming practices
Code
#if (!require("devtools", quietly = TRUE)) {
# install.packages("devtools")
#}
#devtools::install_github("adrientaudiere/MiscMetabar")Code
library(MiscMetabar)
library(targets)
library(metacoder)
library(DT)
library(knitr)
library(ggplot2)
library(patchwork)
library(vegan)
library(FactoMineR)
library(factoextra)
library(DESeq2)
library(indicspecies)
library(sf)
library(adespatial)
library(ade4)
library(adegraphics)
library(kableExtra)
library(formattable)
library(microbiomeMarker)
library(ComplexUpset)Code
d <- tar_read(d_vs_blast)Metadata modification
Code
d@sam_data$rank <- gsub("Weeding", "Grassing", d@sam_data$rank)
d@sam_data$inter_rank <- gsub("Weeding", "Grassing", d@sam_data$inter_rank)
d@sam_data$soil_practice <- gsub("Weeding", "Grassing", d@sam_data$soil_practice)
d@sam_data$terroir <- gsub("Perrier", "Vergèze", d@sam_data$terroir)Spores number
Code
d@sam_data$nb_spores <- rowMeans(cbind(as.numeric(d@sam_data$spores_rep1), as.numeric(d@sam_data$spores_rep2),as.numeric(d@sam_data$spores_rep3)),na.rm=TRUE)Practices
Code
d@sam_data$organic <- ifelse(!d@sam_data$practice %in%
c("Conventional", "Conversion"),
"Organic", "NonOrganic"
)Paired samples (roots, spores from the same samples)
Code
d@sam_data$paired_name <-
apply(as.matrix(d@sam_data[, c(5:16)]), 1, paste,
collapse = "--"
)
d@sam_data$paired_name[d@sam_data$terroir == "Vergèze"] <-
gsub("R", "", d@sam_data$ref_mycea[d@sam_data$terroir == "Vergèze"])Code
tib_sam_data <- as_tibble(d@sam_data) %>%
mutate(across(starts_with("Myc_freq"), as.numeric)) %>%
mutate(across(starts_with("Myc_intensity_colonization"), as.numeric)) %>%
mutate(across(starts_with("Myc_intensity_root"), as.numeric)) %>%
mutate(across(starts_with("Arbuscul_richness"), as.numeric)) %>%
mutate(across(starts_with("Arbuscul_abundance"), as.numeric)) %>%
mutate(across(starts_with("Vesicle_richness"), as.numeric)) %>%
mutate(across(starts_with("Vesicle_abundance"), as.numeric)) %>%
mutate(Myc_freq = rowMeans(select(., starts_with("Myc_freq")), na.rm = TRUE)) %>%
mutate(Myc_intensity_colonization = rowMeans(select(
., starts_with("Myc_intensity_colonization")
), na.rm = TRUE)) %>%
mutate(Myc_intensity_root = rowMeans(select(., starts_with(
"Myc_intensity_root"
)), na.rm = TRUE)) %>%
mutate(Arbuscul_richness = rowMeans(select(., starts_with(
"Arbuscul_richness"
)), na.rm = TRUE)) %>%
mutate(Arbuscul_abundance = rowMeans(select(., starts_with(
"Arbuscul_abundance"
)), na.rm = TRUE)) %>%
mutate(Vesicle_richness = rowMeans(
select(., starts_with(
"Vesicle_richness"
)),
na.rm = TRUE
)) %>%
mutate(Vesicle_abundance = rowMeans(
select(., starts_with(
"Vesicle_abundance"
)),
na.rm = TRUE
)) %>% tibble::column_to_rownames(var = "X")
sample_data(d) <- tib_sam_dataExclude Mycelium
Code
d <- clean_pq(subset_samples(d, compartment != "Mycelium"))
d_vs <- clean_pq(subset_samples(tar_read(d_vs), compartment != "Mycelium"))Summary of final dataset
Code
summary_plot_pq(d) +
labs(title = "Summary of final dataset") +
theme(plot.title = element_text(hjust = 0.5, size = 20, color = "#1e2b4c"))
Summary of bioinformatics pipeline
Global summary
Code
kable(tar_read(track_sequences_samples_clusters))| nb_sequences | nb_clusters | nb_samples | |
|---|---|---|---|
| Paired sequences | 8283716 | 28637 | 194 |
| Paired sequences without chimera | 7777304 | 9791 | 194 |
| Paired sequences without chimera and longer than 200bp | 7736170 | 7522 | 194 |
| ASV denoising | 7533009 | 7522 | 186 |
| OTU after vsearch reclustering at 97% | 7533009 | 1081 | 186 |
| OTU after blast filter without reclustering | 3819912 | 3793 | 182 |
| OTU after blast filter with reclustering | 3822304 | 219 | 182 |
SUBSETING AND FILTERING DATASET
Code
d_muco <- clean_pq(subset_taxa(d, Class_PR2 == "Mucoromycota"))
d_vs_blast <- clean_pq(subset_taxa(tar_read(d_vs_blast), Class_PR2 == "Mucoromycota"))Code
for(i in c(19:46)){
d_muco@sam_data[,i] <- as.numeric(unlist(d_muco@sam_data[,i]))
}Code
d_roots <- clean_pq(subset_samples(d_muco, compartment == "Roots"))
d_spores <- clean_pq(subset_samples(d_muco, compartment == "Spores"))Code
d_rs_merged_samples <- clean_pq(merge_samples2(
d_muco,
"paired_name",
default_fun =
function(x){MiscMetabar::diff_fct_diff_class(x, na.rm = T)}
))Code
d_rs <- clean_pq(subset_samples_pq(
d_rs_merged_samples,
sample_sums(d_rs_merged_samples) > 2000
))TAXONOMICAL ANALYSES
Code
# Figure S8
multitax_bar_pq(d_vs, "Supergroup_PR2",
"Subdivision_PR2", "Class_PR2",
nb_seq = TRUE
) +
ggtitle("Number of sequences (log10) including non-Mucoromycota")
Number of sequences
Code
# Figure S6
dir.create("output/krona")Warning in dir.create("output/krona"): 'output/krona' existe déjàCode
krona(
d_vs,
ranks = c(11:19),
"output/krona/OTU_post_clustered_allEuk",
name = "OTU_post_clustered_allEuk"
)
krona(
d,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjam_AM",
name = "OTU_filtOnMaarjam_AM"
)
krona(
d_muco,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjamAndPr2_AM",
name = "OTU_filtOnMaarjamAndPr2_AM"
)
merge_krona(
c(
"output/krona/OTU_post_clustered_allEuk",
"output/krona/OTU_filtOnMaarjam_AM",
"output/krona/OTU_filtOnMaarjamAndPr2_AM"
),
"output/krona/Euk_to_AM_filtering_nb_seq.html"
)Number of OTU
Code
# Figure S7
krona(
d_vs,
ranks = c(11:19),
"output/krona/OTU_post_clustered_allEuk_Nbotu",
name = "OTU_post_clustered_allEuk",
nb_seq = F
)
krona(
d,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjam_AM_Nbotu",
name = "OTU_filtOnMaarjam_AM",
nb_seq = F
)
krona(
d_muco,
ranks = c(11:19),
"output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu",
name = "OTU_filtOnMaarjamAndPr2_AM",
nb_seq = F
)
merge_krona(
c(
"output/krona/OTU_post_clustered_allEuk_Nbotu",
"output/krona/OTU_filtOnMaarjam_AM_Nbotu",
"output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu"
),
"output/krona/Euk_to_AM_filtering_Nbotu.html"
)NUMBER OF SPORES ANALYSES
Code
psm_otu <- psmelt(as_binary_otu_table(d_spores)) |>
filter(Abundance > 0) |>
group_by(Sample) |>
summarise(
"Abundance" = sum(Abundance),
"region" = region[1],
"nb_spores" = nb_spores[1]
)Code
psm_res <- psmelt_samples_pq(d_spores) %>%
mutate(nb_spores_log10 = log10(nb_spores))Spores count <-> Alpha diversity
Code
# Fig S2
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_spores, rngseed = 221013)) %>% mutate(nb_spores_log10 = log10(nb_spores))
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Abundance, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_0, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_1, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_2, type = "non-parametric")
Spores count <-> Practice
Code
# Figure S5
ggstatsplot::ggbetweenstats(psm_res, practice, nb_spores_log10, type = "non-parametric") 
Spores count <-> Terroir
Code
# Figure S4
ggstatsplot::ggbetweenstats(
mutate(psm_res,
terroir = reorder(psm_res$terroir, psm_res$nb_spores)
),
terroir,
nb_spores_log10,
type = "non-parametric",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist"
) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
Spores count <-> Mycorrhization rate
Code
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_freq, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_intensity_colonization, type = "non-parametric") +
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Arbuscul_abundance, type = "non-parametric")
MYCORRHIZATION MEASURE
Mycorrhization rate <-> Alpha diversity
F = Myc_freq
Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))
# Figure S3a
(ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_2, type = "non-parametrique"))
M = Myc_intensity_colonization
Code
# Figure S3b
(ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_2, type = "non-parametrique"))
A = Arbuscul_abundance
Code
# Figure S3c
(ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_2, type = "non-parametrique")) /
(ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_2, type = "non-parametrique"))
Mycorrhization rate <-> terroir
Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)
Mycorrhization rate <-> practice
Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)
ECOLOGICAL ANALYSES : ALPHA DIVERSITY
Compartment
Code
psm_res <- psmelt_samples_pq(d_muco)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_muco, rngseed = 221013))Code
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_0,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_1,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
compartment,
Hill_2,
type = "non-parametrique"
)
Code
ggplot(psm_res_rarefied, aes(y=Hill_0, x=compartment, color=practice))+
geom_point() +
geom_line(aes(group = paired_name)) +
facet_wrap(~terroir)
Terroir
Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))Code
# Figure 4
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0)
),
terroir,
Hill_0,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1)
),
terroir,
Hill_1,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res_rarefied,
terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2)
),
terroir,
Hill_2,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2))),levels(as.factor(psm_res_rarefied$terroir)))]) + coord_flip()) + plot_layout(axes = "collect")
Code
# Figure S10
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_0)
),
terroir,
Hill_0,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_0))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_1)
),
terroir,
Hill_1,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_1))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) /
(ggstatsplot::ggbetweenstats(
mutate(
psm_res,
terroir = reorder(psm_res$terroir, psm_res$Hill_2)
),
terroir,
Hill_2,
type = "non-parametrique",
centrality.plotting = F,
package = "ggthemes",
palette = "stata_economist",
point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
0.9, size = 3, stroke = 0, na.rm = TRUE),
violin.args = list(width = 0, linewidth = 0)
)+
scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_2))),levels(as.factor(psm_res$terroir)))]) + coord_flip()) + plot_layout(axes = "collect")
Code
# Fig_7a
p_accu_balanced_mod_terroir <-
MiscMetabar::accu_plot_balanced_modality(d_rs,
"terroir",
nperm = 999,
step = 300,
quantile_prob=0.9,
sample.size = 10000,
progress_bar = FALSE)Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3
Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 3Warning in vegan::rarefy(as.matrix(unclass(x[i, ])), n, se = TRUE): most
observed count data have counts 1, but smallest count is 17Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyéCode
p_accu_balanced_mod_terroir +
xlab("Number of sequences") +
ylab("Number of OTUs") +
ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod_terroir$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod_terroir$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,25000)) + theme(legend.position = "none" )Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_line()`).
Practice
Code
# Figure 5
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_0,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 0 (richness)") +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_1,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 1 (~Shannon)") +
ggstatsplot::ggbetweenstats(psm_res_rarefied,
practice,
Hill_2,
centrality.plotting = F,
type = "non-parametrique"
) + ylab("Hill 2 (~Simpson)") + plot_layout(axes = "collect")
Code
# Figure S11
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_0,
centrality.plotting = F,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_1,
centrality.plotting = F,
type = "non-parametrique"
) +
ggstatsplot::ggbetweenstats(psm_res,
practice,
Hill_2,
centrality.plotting = F,
type = "non-parametrique"
)
Code
# Fig_7b
p_accu_balanced_mod <-
MiscMetabar::accu_plot_balanced_modality(d_rs,
"practice",
nperm = 999,
step = 300,
quantile_prob=0.90,
progress_bar = FALSE)Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyéCode
p_accu_balanced_mod +
xlab("Number of sequences") +
ylab("Number of OTUs") +
ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,18000)) + theme(legend.position = "none" )Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_line()`).
ZOOM ON LOW-DIVERSITY SAMPLES
Code
otu_hill <- vegan::renyi(d_rs@otu_table, scale = c(0, 1, 2), hill = TRUE)
low_div_samples <-
rownames(otu_hill[otu_hill$`0` < 10 |
otu_hill$`1` < 2 |
otu_hill$`2` < 2, ])
low_div_samples [1] "2019--Languedoc--Cevennes--St Maurice de Cazevieille--Conventional--Conventional--4.247388--44.035228--Conventional--Herbicide--Scratching--Herbicide"
[2] "PE13"
[3] "2023--Champagne--Côte des Blancs--Chetillons--Chetillons Desh sous rang_plein chimique--Conventional--4.034283--48.945738--Conventional--Herbicide--Herbicide--Herbicide"
[4] "2023--Champagne--Côte des Blancs--Chetillons--Chetillons Ecorces--Conventional--4.034333--48.945724--Conventional--Herbicide--Grassing--Herbicide"
[5] "2023--Champagne--Côte des Blancs--Cramant--BIONNES JEB--Conversion--4.018985--48.988723--Conversion--Scratching--Grassing--Grassing"
[6] "2023--Champagne--Côte des Blancs--Cramant--BATEAU AN39 JEB--Conversion--4.021354--48.995715--Conversion--Scratching--Grassing--Grassing"
[7] "2023--Champagne--Côte des Blancs--Cramant--TDB AM29 JEB--Conversion--4.016860--48.996831--Conversion--Scratching--Grassing--Grassing"
[8] "2023--Champagne--Côte des Blancs--Mireaux--JEB VOISIN 1 MIREAUX DESSOUS--Conventional--3.994335--48.984188--Conventional--Herbicide--Herbicide--Herbicide"
[9] "2023--Champagne--Côte des Blancs--Avize--AVIZE Desh chimique en plein--Conventional--4.014914--48.976412--Conventional--Herbicide--Herbicide--Herbicide"
[10] "2023--Champagne--Côte des Blancs--Avize--AVIZE Travail meca sous rang--Organic--4.014914--48.976351--Organic grassed--Scratching--Grassing--Grassing"
[11] "2020--Champagne--Côte des Bar--Les Riceys--Conventional--Conventional--4.392132--48.009899--Conventional--Herbicide--Scratching--Herbicide"
[12] "2020--Bordelais--Langoiran--Langoiran--Conventional chimique--Conventional---0.356503--44.713204--Conventional--Herbicide--Herbicide--Herbicide"
[13] "2020--Champagne--Vallee de la Marne--Vincelles--Conventional--Conventional--3.931795--49.077467--Conventional--Herbicide--Scratching--Herbicide"
[14] "2020--Camargue--Montcalm--Montcalm--Conventional--Conventional--4.320967--43.584075--Conventional--Scratching--Scratching--Scratching" Code
sum(otu_hill$`0` < 10)[1] 11Code
sum(otu_hill$`1` < 2)[1] 8Code
sum(otu_hill$`2` < 2)[1] 11Code
length(low_div_samples)[1] 14Code
low_div_samples_table <- tibble(cbind(as_tibble(d_rs@sam_data[low_div_samples, c("terroir", "practice", "rank", "inter_rank", "compartment")]), apply(otu_hill[low_div_samples, ], 2, round, digits = 2), sample_sums(d_rs)[low_div_samples]))Warning in class(x) <- tibble_class: La class(x) est une chaîne multiple
("tbl_df", "tbl", ...) ; le résultat ne sera plus un objet S4Code
low_div_samples_table$compartment <- ifelse(is.na(low_div_samples_table$compartment), "Spores + Roots", "Spores")
colnames(low_div_samples_table) <- c("Terroir", "Global practice", "Rank practice", "Inter rank practice", "Compartment", "Hill 0 (richness)", "Hill 1 (~Shannon)", "Hill 2 (~Simpson)", "nb_seq")
low_div_samples_table <- arrange(low_div_samples_table, as.numeric(`Hill 0 (richness)`))
low_div_samples_table <-
rbind(
low_div_samples_table,
c(
"",
"",
"",
"MEAN",
"(low -diversity samples)",
round(mean(low_div_samples_table$`Hill 0 (richness)`), 2),
round(mean(low_div_samples_table$`Hill 1 (~Shannon)`), 2),
round(mean(low_div_samples_table$`Hill 2 (~Simpson)`), 2),
round(mean(low_div_samples_table$nb_seq), 2)
),
c(
"",
"",
"",
"MEAN",
"(all samples)",
round(mean(otu_hill$`0`), 2),
round(mean(otu_hill$`1`), 2),
round(mean(otu_hill$`2`), 2),
round(mean(sample_sums(d_rs)), 2)
),
c(
"",
"",
"",
"MAX",
"(all samples)",
round(max(otu_hill$`0`), 2),
round(max(otu_hill$`1`), 2),
round(max(otu_hill$`2`), 2),
round(max(sample_sums(d_rs)), 2)
)
)
# Table 3
kbl(low_div_samples_table) |>
kable_classic(full_width = F, html_font = "Cambria")| Terroir | Global practice | Rank practice | Inter rank practice | Compartment | Hill 0 (richness) | Hill 1 (~Shannon) | Hill 2 (~Simpson) | nb_seq |
|---|---|---|---|---|---|---|---|---|
| Côte des Blancs | Conversion | Scratching | Grassing | Spores + Roots | 3 | 1.87 | 1.68 | 42775 |
| Côte des Blancs | Conventional | Herbicide | Grassing | Spores + Roots | 4 | 1.71 | 1.39 | 28481 |
| Langoiran | Conventional | Herbicide | Herbicide | Spores + Roots | 4 | 2.38 | 1.87 | 12643 |
| Côte des Blancs | Organic | Scratching | Grassing | Spores + Roots | 5 | 1.99 | 1.49 | 35728 |
| Vergèze | Organic | Scratching | Ploughing | Spores | 6 | 4.33 | 3.65 | 6768 |
| Côte des Blancs | Conventional | Herbicide | Herbicide | Spores + Roots | 6 | 1.34 | 1.12 | 41991 |
| Côte des Blancs | Conversion | Scratching | Grassing | Spores + Roots | 7 | 2.74 | 2.13 | 49313 |
| Côte des Blancs | Conventional | Herbicide | Herbicide | Spores + Roots | 8 | 1.68 | 1.26 | 46114 |
| Côte des Blancs | Conversion | Scratching | Grassing | Spores + Roots | 8 | 1.17 | 1.05 | 45845 |
| Vallee de la Marne | Conventional | Herbicide | Scratching | Spores + Roots | 9 | 2.7 | 1.78 | 11777 |
| Montcalm | Conventional | Scratching | Scratching | Spores | 9 | 3.45 | 2.94 | 24675 |
| Côte des Blancs | Conventional | Herbicide | Herbicide | Spores + Roots | 10 | 1.45 | 1.16 | 45017 |
| Cevennes | Conventional | Herbicide | Scratching | Spores + Roots | 12 | 1.91 | 1.35 | 93267 |
| Côte des Bar | Conventional | Herbicide | Scratching | Spores + Roots | 17 | 2.43 | 1.51 | 28971 |
| MEAN | (low -diversity samples) | 7.71 | 2.22 | 1.74 | 36668.93 | |||
| MEAN | (all samples) | 19.87 | 6.23 | 4.35 | 48209.53 | |||
| MAX | (all samples) | 47 | 15.06 | 9.04 | 101104 |
ECOLOGICAL ANALYSES : BETA DIVERSITY
Compartment
Code
# Figure S9
p <- phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), "matrix"),
method = "robust.aitchison"
)),
color = "region",
shape = "compartment"
) +
geom_point(size = 3) +
stat_ellipse(inherit.aes = F, aes(x = DCA1, y = DCA2, linetype = compartment))Warning in phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), : `Ordination species/OTU/taxa
coordinate indices did not match `physeq` index names. Setting corresponding
coordinates to NULL.Code
p + geom_line(data = p$data, aes(group = paired_name))
Practice and terroir
Code
p <-
phyloseq::plot_ordination(d_rs,
ordinate(d_rs,
method = "NMDS",
distance = "bray"
),
color = "terroir",
shape = "practice"
) + facet_wrap("region") + scale_shape_manual(values=c(25,23,24))Square root transformation
Wisconsin double standardization
Run 0 stress 0.2485558
Run 1 stress 0.2596844
Run 2 stress 0.2663072
Run 3 stress 0.2567211
Run 4 stress 0.256791
Run 5 stress 0.2626177
Run 6 stress 0.2636841
Run 7 stress 0.2612171
Run 8 stress 0.2652946
Run 9 stress 0.262289
Run 10 stress 0.2523773
Run 11 stress 0.2583239
Run 12 stress 0.2599351
Run 13 stress 0.2562442
Run 14 stress 0.409148
Run 15 stress 0.2593146
Run 16 stress 0.2658447
Run 17 stress 0.2580753
Run 18 stress 0.2626311
Run 19 stress 0.261243
Run 20 stress 0.2584539
*** Best solution was not repeated -- monoMDS stopping criteria:
2: no. of iterations >= maxit
18: stress ratio > sratmaxCode
p_nmds <- p + geom_point(
data = select(p$data, -c(region)),
fill = "grey40",
color = "grey20",
size = 1,
alpha = 0.5
) +
geom_point(size = 3, aes(fill=terroir)) +
ggtitle("NMDS on bray distance") +
scale_fill_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
"#e15968","#784b43","#c200ab","#00b3f0","#953726",
"#e09199","#379d30","#ed5019","#ff63b0")) +
scale_color_manual(values=rep("grey20",14))
# Figure 9
p_nmds
Code
# Figure 6a
upset_pq(d_rs, "practice", taxa_fill = "Family",
set_sizes=(
upset_set_size()
+ geom_text(aes(label=after_stat(count)), hjust=-0.3, color="white", stat='count')
))
Code
ggvenn_pq(d_rs, "practice")
Code
# Figure 6b
upset_pq(d_rs, "terroir", taxa_fill = "Family", min_size = 2, height_ratio = 0.6
,set_sizes=(
upset_set_size()
+ geom_label(aes(label=after_stat(count)), hjust=-0.3, size=2.5, stat='count')
))Warning: Removed 61 rows containing non-finite outside the scale range
(`stat_count()`).
Global beta-div analysis
Spatial data
Code
library(geodist)
lat_lon <- as_tibble(d_rs@sam_data[, c("lat", "long")])
lat_lon$lat <- as.numeric(lat_lon$lat)
lat_lon$long <- as.numeric(lat_lon$long)
dist_spatial_meter <- as.dist(geodist(lat_lon, measure = "geodesic"),
upper = FALSE
)Code
lon_lat_rs <- d_rs@sam_data[, c("long", "lat")]
lon_lat_rs$long <- as.numeric(lon_lat_rs$long)
lon_lat_rs$lat <- as.numeric(lon_lat_rs$lat)
MEM <- dbmem(lon_lat_rs, MEM.autocor = "non-null")Code
test_MEM <- moran.randtest(MEM, nrepet = 999)
test_MEM$pvalue_adjust <- p.adjust(test_MEM$pvalue, method = "BH")Code
d_rs@sam_data$MEM_1 <- MEM[, 1]
d_rs@sam_data$MEM_2 <- MEM[, 2]
res_ado_spatial_robAit <- adonis_pq(
d_rs,
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "robust.aitchison"
)
res_ado_spatial_robAit_rarefy <- adonis_pq(
rarefy_even_depth(d_rs, rngseed = 626),
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "robust.aitchison"
)
res_ado_spatial_bray <- adonis_pq(
d_rs,
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "bray"
)
res_ado_spatial_bray_rarefy <- adonis_pq(
rarefy_even_depth(d_rs, rngseed = 626),
"MEM_1 + MEM_2 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "bray"
)
df <- data.frame(res_ado_spatial_bray)
df$names <-
factor(
rownames(df),
levels = c(
"sample_size",
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p1 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on Bray")
df <- data.frame(res_ado_spatial_robAit)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p2 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on robust Aitchison")
df <- data.frame(res_ado_spatial_bray_rarefy)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p3 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on Bray after rarefaction")
df <- data.frame(res_ado_spatial_robAit_rarefy)
df$names <-
factor(
rownames(df),
levels = c(
"MEM_1",
"MEM_2",
"practice",
"inter_rank",
"rank",
"terroir",
"Residual",
"Total"
)
)
df <- df %>% filter(!names == "Total")
p4 <-
ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
geom_bar(stat = "identity") +
scale_fill_viridis_c() +
geom_text(
label = case_when(
df$Pr..F. < 0.001 ~ "***",
df$Pr..F. < 0.01 ~ "**",
df$Pr..F. < 0.05 ~ "*",
df$Pr..F. > 0.05 ~ "ns"
),
nudge_y = 0.01
) +
ggtitle("Permanova on robust Aitchison after rarefaction")
(p1 + ylim(0, 0.12) + p2 + ylim(0, 0.12)) /
(p3 + ylim(0, 0.12) + p4 + ylim(0, 0.12))Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.16217 0.012475 0.6941 0.7615
Residuals 61 1.09633 0.017973 Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.05732 0.028658 2.1729 0.1213
Residuals 72 0.94962 0.013189 Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.16217 0.012475 0.6941 0.7615
Residuals 61 1.09633 0.017973 Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.05732 0.028658 2.1729 0.1213
Residuals 72 0.94962 0.013189 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 393.8 30.2920 9.5001 2.306e-10 ***
Residuals 61 194.5 3.1886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 39.37 19.685 1.4366 0.2445
Residuals 72 986.60 13.703 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 393.8 30.2920 9.5001 2.306e-10 ***
Residuals 61 194.5 3.1886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 39.37 19.685 1.4366 0.2445
Residuals 72 986.60 13.703 Soil properties
Code
soil_prop <- as_tibble(d_rs@sam_data) |>
select(paired_name, practice, organic, terroir, Coarse_sand:CEC)
soil_prop_num <- soil_prop |>
select(-all_of(c("practice", "organic", "terroir", "CaO", "Total_sand", "Total_filt"))) |>
tibble::column_to_rownames("paired_name") |>
mutate(across(everything(), as.numeric)) |>
tidyr::drop_na()Code
pca_test_res <- PCAtest::PCAtest(soil_prop_num, varcorr = T, plot = F, counter = F)
Sampling bootstrap replicates... Please wait
Calculating confidence intervals of empirical statistics... Please wait
Sampling random permutations... Please wait
Comparing empirical statistics with their null distributions... Please wait
========================================================
Test of PCA significance: 17 variables, 44 observations
1000 bootstrap replicates, 1000 random permutations
========================================================
Empirical Psi = 55.8085, Max null Psi = 9.0868, Min null Psi = 3.6595, p-value = 0
Empirical Phi = 0.4530, Max null Phi = 0.1828, Min null Phi = 0.1160, p-value = 0
Empirical eigenvalue #1 = 7.28487, Max null eigenvalue = 2.88796, p-value = 0
Empirical eigenvalue #2 = 3.15651, Max null eigenvalue = 2.39696, p-value = 0
Empirical eigenvalue #3 = 2.59627, Max null eigenvalue = 2.09788, p-value = 0
Empirical eigenvalue #4 = 1.16046, Max null eigenvalue = 1.87236, p-value = 1
Empirical eigenvalue #5 = 0.93253, Max null eigenvalue = 1.59515, p-value = 1
Empirical eigenvalue #6 = 0.74202, Max null eigenvalue = 1.48039, p-value = 1
Empirical eigenvalue #7 = 0.39204, Max null eigenvalue = 1.32508, p-value = 1
Empirical eigenvalue #8 = 0.24218, Max null eigenvalue = 1.18435, p-value = 1
Empirical eigenvalue #9 = 0.19075, Max null eigenvalue = 1.0485, p-value = 1
Empirical eigenvalue #10 = 0.08962, Max null eigenvalue = 0.99123, p-value = 1
Empirical eigenvalue #11 = 0.07595, Max null eigenvalue = 0.85026, p-value = 1
Empirical eigenvalue #12 = 0.05693, Max null eigenvalue = 0.7427, p-value = 1
Empirical eigenvalue #13 = 0.0369, Max null eigenvalue = 0.65767, p-value = 1
Empirical eigenvalue #14 = 0.02507, Max null eigenvalue = 0.55583, p-value = 1
Empirical eigenvalue #15 = 0.01767, Max null eigenvalue = 0.48237, p-value = 1
Empirical eigenvalue #16 = 0.00018, Max null eigenvalue = 0.40241, p-value = 1
Empirical eigenvalue #17 = 5e-05, Max null eigenvalue = 0.35079, p-value = 1
PC 1 is significant and accounts for 42.9% (95%-CI:38.8-49.7) of the total variation
PC 2 is significant and accounts for 18.6% (95%-CI:16.6-25) of the total variation
PC 3 is significant and accounts for 15.3% (95%-CI:10.5-17.9) of the total variation
The first 3 PC axes are significant and account for 76.7% of the total variation
Variables 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, and 17 have significant loadings on PC 1
Variables 4, 8, 10, 12, and 14 have significant loadings on PC 2
Variables 5, 6, and 7 have significant loadings on PC 3
Variables 8, 10, 13, 15, and 17 have significant correlations with PC 1
Variables 8, 10, 12, and 14 have significant correlations with PC 2
Variables , and have significant correlations with PC 3Code
pval <- c()
for (i in seq_len(length(pca_test_res$`Percentage of variation of empirical PC's`))) {
obs <- pca_test_res$`Percentage of variation of empirical PC's`[i]
null_model <- pca_test_res$`Percentage of variation of randomized data`[, i]
pval[i] <- (sum(obs < null_model) + 1) / (1 + nrow(pca_test_res$`Percentage of variation of randomized data`))
}
pca_test_pval_adj <- p.adjust(pval, method = "BH")
pca_test_pval_adj [1] 0.005661006 0.005661006 0.005661006 1.000000000 1.000000000 1.000000000
[7] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
[13] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000Code
res_pca <- dudi.pca(soil_prop_num, scannf = F, nf = 4)Code
p_pca_variable <- fviz_pca_var(res_pca,
col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
p_pca_variable_axe3_4 <- fviz_pca_var(res_pca,
col.var = "cos2", axes = c(3, 4),
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
res_pca_var <- get_pca_var(res_pca)
# Fig. 2
(p_pca_variable + p_pca_variable_axe3_4) / free((factoextra::fviz_eig(res_pca)) +
(ggcorrplot::ggcorrplot(res_pca_var$cor)))
Code
practice <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
x[5]
}))
terroir <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
x[2]
}))
p_pca_ind_terroir <- fviz_pca_ind(res_pca,
habillage = terroir,
label = FALSE,
repel = TRUE,
mean.point = FALSE,
addEllipses = TRUE,
pointsize = 3,
ellipse.type = "convex"
) +
scale_shape_manual(values = rep(16, length(unique(terroir))))
p_pca_ind_terroir$data$terroir <- terroir
p_pca_ind_terroir$data$practice <- practice
p_pca_ind_terroir_practice <- p_pca_ind_terroir + geom_point(aes(x = x, y = y),
shape = case_when(
practice == "Organic" ~ 2,
practice == "Conventional" ~ 6,
practice == "Conversion" ~ 5
), size = 3.5,
color = rgb(0, 0, 0, 0.7)
)
# value
kw1_terroir <- kruskal.test(res_pca$tab[, 1], terroir)
kw2_terroir <- kruskal.test(res_pca$tab[, 2], terroir)
kw3_terroir <- kruskal.test(res_pca$tab[, 3], terroir)
kw1_practice <- kruskal.test(res_pca$tab[, 1], practice)
kw2_practice <- kruskal.test(res_pca$tab[, 2], practice)
kw3_practice <- kruskal.test(res_pca$tab[, 3], practice)
# Figure 3
p_pca_ind_terroir_practice +
annotate("text", x = 2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - terroir): pval=", round(as.numeric(kw1_terroir$p.value),4)), size = 3) +
annotate("text", x = 2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - terroir): pval=", round(as.numeric(kw2_terroir$p.value),4)), size = 3) +
annotate("text", x = 2, y = 4, label = paste("Kruskal-Wallis (axe3 - terroir): pval=", round(as.numeric(kw3_terroir$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - practice): pval=", round(as.numeric(kw1_practice$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - practice): pval=", round(as.numeric(kw2_practice$p.value),4)), size = 3) +
annotate("text", x = -2, y = 4, label = paste("Kruskal-Wallis (axe3 - practice): pval=", round(as.numeric(kw3_practice$p.value),4)), size = 3)Warning: Removed 21 rows containing missing values or values outside the scale range
(`geom_point()`).
Code
# Figure S1
ggplot(d_rs@sam_data, aes(x = practice, y = as.numeric(Cu), fill = practice)) +
geom_violin() +
geom_jitter()Warning: Removed 23 rows containing non-finite outside the scale range
(`stat_ydensity()`).Warning: Removed 23 rows containing missing values or values outside the scale range
(`geom_point()`).
Code
kruskal.test(as.numeric(d_rs@sam_data$Cu), d_rs@sam_data$practice)
Kruskal-Wallis rank sum test
data: as.numeric(d_rs@sam_data$Cu) and d_rs@sam_data$practice
Kruskal-Wallis chi-squared = 0.72149, df = 2, p-value = 0.6972Code
d_with_pca <-
clean_pq(subset_samples_pq(d_rs, d_rs@sam_data$paired_name %in% rownames(soil_prop_num)))
res_pca_ind <- get_pca_ind(res_pca)
if (!identical(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name), sort(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name)))) {
stop("ERROR")
}
d_with_pca <- add_info_to_sam_data(d_with_pca, res_pca_ind$coord)Code
dim_df <- as_tibble(d_with_pca@sam_data) |>
select(c(starts_with("Dim"), "nb_seq", "nb_otu"))
cor_results <- correlation::correlation(dim_df)
cor_results %>%
summary(redundant = TRUE) %>%
plot()
Code
lon_lat_with_pca <- d_with_pca@sam_data[, c("long", "lat")]
lon_lat_with_pca$long <- as.numeric(lon_lat_with_pca$long)
lon_lat_with_pca$lat <- as.numeric(lon_lat_with_pca$lat)
MEM_with_pca <- dbmem(lon_lat_with_pca, MEM.autocor = "non-null")
d_with_pca@sam_data$MEM_1 <- MEM_with_pca[, 1]
d_with_pca@sam_data$MEM_2 <- MEM_with_pca[, 2]Code
res_ado_spatial_soil_bray <-
adonis_pq(
d_with_pca,
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE
)
res_ado_spatial_soil_robAit <-
adonis_pq(
d_with_pca,
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = TRUE,
dist_method = "robust.aitchison"
)
res_ado_spatial_soil_bray_rarefy <-
adonis_pq(
rarefy_even_depth(d_with_pca, rngseed = 626),
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE
)
res_ado_spatial_soil_robAit_rarefy <-
adonis_pq(
rarefy_even_depth(d_with_pca, rngseed = 626),
"MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
correction_for_sample_size = FALSE,
dist_method = "robust.aitchison"
)
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.05732 0.028658 2.1729 0.1213
Residuals 72 0.94962 0.013189 Code
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 0.16217 0.012475 0.6941 0.7615
Residuals 61 1.09633 0.017973 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 39.37 19.685 1.4366 0.2445
Residuals 72 986.60 13.703 Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 13 393.8 30.2920 9.5001 2.306e-10 ***
Residuals 61 194.5 3.1886
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
# TABLE 2
kbl(res_ado_spatial_bray) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| sample_size | 1 | 1.0052204 | 0.0670420 | 7.319743 | 0.001 |
| MEM_1 | 1 | 1.5103278 | 0.1007296 | 10.997797 | 0.001 |
| MEM_2 | 1 | 0.9164972 | 0.0611247 | 6.673684 | 0.001 |
| practice | 2 | 0.4789872 | 0.0319455 | 1.743928 | 0.036 |
| inter_rank | 3 | 0.5442769 | 0.0362999 | 1.321092 | 0.128 |
| rank | 2 | 0.6046954 | 0.0403295 | 2.201614 | 0.004 |
| terroir | 13 | 2.9300449 | 0.1954160 | 1.641215 | 0.001 |
| Residual | 51 | 7.0038313 | 0.4671126 | NA | NA |
| Total | 74 | 14.9938812 | 1.0000000 | NA | NA |
Code
kbl(res_ado_spatial_soil_bray) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| sample_size | 1 | 0.9095040 | 0.1022131 | 6.3231283 | 0.001 |
| MEM_1 | 1 | 0.2091359 | 0.0235034 | 1.4539720 | 0.149 |
| MEM_2 | 1 | 1.0896243 | 0.1224557 | 7.5753755 | 0.001 |
| Dim.1 | 1 | 0.3845781 | 0.0432202 | 2.6736956 | 0.007 |
| Dim.2 | 1 | 0.2413033 | 0.0271185 | 1.6776088 | 0.080 |
| Dim.3 | 1 | 0.0934249 | 0.0104994 | 0.6495160 | 0.783 |
| practice | 2 | 0.2815596 | 0.0316426 | 0.9787409 | 0.508 |
| inter_rank | 2 | 0.1453520 | 0.0163352 | 0.5052642 | 0.957 |
| rank | 2 | 0.4576684 | 0.0514343 | 1.5909201 | 0.070 |
| terroir | 11 | 2.2092074 | 0.2482782 | 1.3962758 | 0.022 |
| Residual | 20 | 2.8767531 | 0.3232993 | NA | NA |
| Total | 43 | 8.8981110 | 1.0000000 | NA | NA |
Code
# TABLE S2
kbl(res_ado_spatial_bray_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 1.5224735 | 0.1079197 | 11.273304 | 0.001 |
| MEM_2 | 1 | 1.1172045 | 0.0791924 | 8.272450 | 0.001 |
| practice | 2 | 0.5880731 | 0.0416852 | 2.177222 | 0.011 |
| inter_rank | 3 | 0.3890128 | 0.0275750 | 0.960161 | 0.478 |
| rank | 2 | 0.5861573 | 0.0415495 | 2.170130 | 0.016 |
| terroir | 13 | 2.8818779 | 0.2042804 | 1.641473 | 0.003 |
| Residual | 52 | 7.0226637 | 0.4977978 | NA | NA |
| Total | 74 | 14.1074628 | 1.0000000 | NA | NA |
Code
# TABLE S3
kbl(res_ado_spatial_soil_bray_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 0.3924837 | 0.0479272 | 2.8426837 | 0.014 |
| MEM_2 | 1 | 1.0983511 | 0.1341225 | 7.9551458 | 0.001 |
| Dim.1 | 1 | 0.4058745 | 0.0495624 | 2.9396712 | 0.008 |
| Dim.2 | 1 | 0.2866462 | 0.0350031 | 2.0761235 | 0.054 |
| Dim.3 | 1 | 0.1150841 | 0.0140532 | 0.8335323 | 0.581 |
| practice | 2 | 0.2932271 | 0.0358067 | 1.0618937 | 0.388 |
| inter_rank | 2 | 0.1439472 | 0.0175778 | 0.5212911 | 0.945 |
| rank | 2 | 0.4694158 | 0.0573216 | 1.6999441 | 0.061 |
| terroir | 11 | 2.0847081 | 0.2545690 | 1.3726490 | 0.056 |
| Residual | 21 | 2.8994281 | 0.3540566 | NA | NA |
| Total | 43 | 8.1891660 | 1.0000000 | NA | NA |
Code
# TABLE S4
kbl(res_ado_spatial_robAit_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 118.30117 | 0.0313656 | 2.7789327 | 0.001 |
| MEM_2 | 1 | 110.52944 | 0.0293050 | 2.5963722 | 0.002 |
| practice | 2 | 117.62397 | 0.0311860 | 1.3815125 | 0.043 |
| inter_rank | 3 | 93.00477 | 0.0246586 | 0.7282373 | 0.965 |
| rank | 2 | 154.47993 | 0.0409577 | 1.8143917 | 0.002 |
| terroir | 13 | 964.07334 | 0.2556078 | 1.7420303 | 0.001 |
| Residual | 52 | 2213.67750 | 0.5869192 | NA | NA |
| Total | 74 | 3771.69011 | 1.0000000 | NA | NA |
Code
# TABLE S5
kbl(res_ado_spatial_soil_robAit_rarefy) |>
kable_classic(full_width = F, html_font = "Cambria")| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| MEM_1 | 1 | 112.72637 | 0.0339610 | 1.6387773 | 0.008 |
| MEM_2 | 1 | 112.45902 | 0.0338805 | 1.6348907 | 0.006 |
| Dim.1 | 1 | 191.81334 | 0.0577875 | 2.7885165 | 0.001 |
| Dim.2 | 1 | 93.30719 | 0.0281106 | 1.3564679 | 0.087 |
| Dim.3 | 1 | 112.44158 | 0.0338752 | 1.6346370 | 0.014 |
| practice | 2 | 128.09782 | 0.0385919 | 0.9311211 | 0.678 |
| inter_rank | 2 | 98.34228 | 0.0296275 | 0.7148331 | 0.956 |
| rank | 2 | 158.00090 | 0.0476008 | 1.1484814 | 0.175 |
| terroir | 11 | 867.57636 | 0.2613741 | 1.1465934 | 0.057 |
| Residual | 21 | 1444.52439 | 0.4351909 | NA | NA |
| Total | 43 | 3319.28925 | 1.0000000 | NA | NA |
Variance partitionning
Code
# Figure 8a
res_varpart_rarefy <- var_par_rarperm_pq(
physeq = d_with_pca,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Soil" = c("Dim.1", "Dim.2", "Dim.3"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Code
# Figure 8b
res_varpart_rarefy_wo_soil <- var_par_rarperm_pq(
physeq = d_rs,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Code
res_varpart_robAit_rarefy <- var_par_rarperm_pq(
physeq = d_with_pca,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Soil" = c("Dim.1", "Dim.2", "Dim.3"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
dist_method = "robust.aitchison",
nperm = 99,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Code
# Figure S12b
res_varpart_robAit_rarefy_wo_soil <- var_par_rarperm_pq(
physeq = d_rs,
list_component = list(
"Spatial" = c("MEM_1", "MEM_2"),
"Practice" = c("practice", "inter_rank", "rank"),
"Terroir" = c("terroir")
),
dist_method = "robust.aitchison",
nperm = 99,
dbrda_computation = TRUE,
progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)
Differential abundance analysis
Indicspecies
Terroir
Code
# Figure 10
res_mpt_terroir <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$terroir,
control = how(nperm = 9999),
max.order = 3
)
res_mpt_terroir_df <- res_mpt_terroir$sign
res_mpt_terroir_df$p.adj <- p.adjust(res_mpt_terroir_df$p.value, method = "BH")
res_mpt_terroir_df$OTU_names <- rownames(res_mpt_terroir_df)
res_mpt_terroir_df_signif <-
res_mpt_terroir_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
tax_tab <- as.data.frame(d_rs@tax_table)
tax_tab$otu <- rownames(tax_tab)
res_mpt_terroir_df_signif_taxo <- left_join(res_mpt_terroir_df_signif,tax_tab, by = join_by("OTU_names" == "otu"))
ggplot(
res_mpt_terroir_df_signif_taxo,
aes(
x = OTU_names,
y = name,
alpha = value,
color = Genus
)
) +
geom_point(size = 6) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))+
scale_alpha(range = c(0, 1))
Code
res_mpt_terroir_rg <-
multipatt(as.matrix(rarefy_even_depth(d_rs, rngseed = 626)@otu_table),
d_rs@sam_data$terroir,
control = how(nperm = 9999),
max.order = 3,
func = "r.g"
)
res_mpt_terroir_rg_df <- res_mpt_terroir_rg$sign
res_mpt_terroir_rg_df$p.adj <- p.adjust(res_mpt_terroir_rg_df$p.value, method = "BH")
res_mpt_terroir_rg_df$OTU_names <- rownames(res_mpt_terroir_rg_df)
res_mpt_terroir_rg_df_signif <-
res_mpt_terroir_rg_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
ggplot(
res_mpt_terroir_rg_df_signif,
aes(
x = OTU_names,
y = name,
alpha = value,
color = stat
)
) +
geom_point(size = 4) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
))
practice
Code
res_mpt_practice <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$practice,
control = how(nperm = 9999),
max.order = 3
)
res_mpt_practice_df <- res_mpt_practice$sign
res_mpt_practice_df$p.adj <- p.adjust(res_mpt_practice_df$p.value, method = "BH")
res_mpt_practice_df$OTU_names <- rownames(res_mpt_practice_df)
res_mpt_practice_df_signif <-
res_mpt_practice_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_df_signif# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
# OTU_names <chr>, name <chr>, value <dbl>Code
res_mpt_practice_rg <-
multipatt(as.matrix(d_rs@otu_table),
d_rs@sam_data$practice,
control = how(nperm = 9999),
max.order = 3,
func = "r.g"
)
res_mpt_practice_rg_df <- res_mpt_practice_rg$sign
res_mpt_practice_rg_df$p.adj <- p.adjust(res_mpt_practice_rg_df$p.value, method = "BH")
res_mpt_practice_rg_df$OTU_names <- rownames(res_mpt_practice_rg_df)
res_mpt_practice_rg_df_signif <-
res_mpt_practice_rg_df %>%
filter(p.adj < 0.05) %>%
tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_rg_df_signif# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
# OTU_names <chr>, name <chr>, value <dbl>Session information
Code
sessionInfo()R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 11 (bullseye)
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
[5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
[7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Paris
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] geodist_0.0.8 ComplexUpset_1.3.3
[3] microbiomeMarker_1.8.0 formattable_0.2.1
[5] kableExtra_1.4.0 adegraphics_1.0-21
[7] ade4_1.7-22 adespatial_0.3-23
[9] sf_1.0-16 indicspecies_1.7.14
[11] DESeq2_1.42.1 SummarizedExperiment_1.32.0
[13] Biobase_2.62.0 MatrixGenerics_1.14.0
[15] matrixStats_1.3.0 GenomicRanges_1.54.1
[17] GenomeInfoDb_1.38.8 IRanges_2.36.0
[19] S4Vectors_0.40.2 BiocGenerics_0.48.1
[21] factoextra_1.0.7 FactoMineR_2.11
[23] vegan_2.6-6.1 lattice_0.22-6
[25] permute_0.9-7 patchwork_1.2.0
[27] knitr_1.46 DT_0.33
[29] metacoder_0.3.7 targets_1.7.0
[31] MiscMetabar_0.9.2 purrr_1.0.2
[33] dplyr_1.1.4 dada2_1.30.0
[35] Rcpp_1.0.12 ggplot2_3.5.1
[37] phyloseq_1.46.0
loaded via a namespace (and not attached):
[1] igraph_2.0.3 mia_1.10.0
[3] Formula_1.2-5 scater_1.30.1
[5] rematch2_2.1.2 zlibbioc_1.48.2
[7] tidyselect_1.2.1 bit_4.0.5
[9] doParallel_1.0.17 BWStest_0.2.3
[11] clue_0.3-65 rjson_0.2.21
[13] blob_1.2.4 stringr_1.5.1
[15] rngtools_1.5.2 S4Arrays_1.2.1
[17] parallel_4.3.3 RNeXML_2.4.11
[19] correlation_0.8.4 png_0.1-8
[21] cli_3.6.3 ggplotify_0.1.2
[23] CVXR_1.0-12 bayestestR_0.13.2
[25] multtest_2.58.0 MultiAssayExperiment_1.28.0
[27] bluster_1.12.0 BiocNeighbors_1.20.2
[29] TreeSummarizedExperiment_2.10.0 adegenet_2.1.10
[31] mime_0.12 evaluate_0.23
[33] tidytree_0.4.6 ComplexHeatmap_2.18.0
[35] ggh4x_0.2.8 stringi_1.8.4
[37] backports_1.4.1 PMCMRplus_1.9.10
[39] lmerTest_3.1-3 gsl_2.1-8
[41] XML_3.99-0.16.1 Exact_3.2
[43] ggVennDiagram_1.5.2 httpuv_1.6.15
[45] Wrench_1.20.0 paletteer_1.6.0
[47] magrittr_2.0.3 leaps_3.1
[49] splines_4.3.3 jpeg_0.1-10
[51] doRNG_1.8.6 wk_0.9.1
[53] rootSolve_1.8.2.4 ggbeeswarm_0.7.2
[55] DBI_1.2.2 withr_3.0.0
[57] class_7.3-22 systemfonts_1.0.6
[59] htmlwidgets_1.6.4 fs_1.6.4
[61] SuppDists_1.1-9.7 SingleCellExperiment_1.24.0
[63] ggrepel_0.9.5 labeling_0.4.3
[65] SparseArray_1.2.4 cellranger_1.1.0
[67] flashClust_1.01-2 rncl_0.8.7
[69] see_0.8.4 lmom_3.0
[71] effectsize_0.8.8 zoo_1.8-12
[73] XVector_0.42.0 decontam_1.22.0
[75] secretbase_0.5.0 foreach_1.5.2
[77] fansi_1.0.6 caTools_1.18.2
[79] grid_4.3.3 data.table_1.15.4
[81] ggtree_3.10.1 rhdf5_2.46.1
[83] irlba_2.3.5.1 gridGraphics_0.5-1
[85] DescTools_0.99.54 base64url_1.4
[87] lazyeval_0.2.2 yaml_2.3.8
[89] survival_3.6-4 crayon_1.5.3
[91] RColorBrewer_1.1-3 tidyr_1.3.1
[93] later_1.3.2 codetools_0.2-19
[95] base64enc_0.1-3 GlobalOptions_0.1.2
[97] shape_1.4.6.1 limma_3.58.1
[99] estimability_1.5.1 Rsamtools_2.18.0
[101] ggside_0.3.1 foreign_0.8-86
[103] pkgconfig_2.0.3 ggpubr_0.6.0
[105] xml2_1.3.6 GenomicAlignments_1.38.2
[107] aplot_0.2.2 scatterplot3d_0.3-44
[109] ape_5.8 viridisLite_0.4.2
[111] xtable_1.8-4 interp_1.1-6
[113] car_3.1-2 highr_0.10
[115] hwriter_1.3.2.1 plyr_1.8.9
[117] httr_1.4.7 rbibutils_2.2.16
[119] tools_4.3.3 broom_1.0.5
[121] beeswarm_0.4.0 htmlTable_2.4.2
[123] checkmate_2.3.1 nlme_3.1-164
[125] PCAtest_0.0.1 lme4_1.1-35.3
[127] digest_0.6.36 numDeriv_2016.8-1.1
[129] adephylo_1.1-16 Matrix_1.6-5
[131] farver_2.1.2 reshape2_1.4.4
[133] yulab.utils_0.1.4 viridis_0.6.5
[135] DirichletMultinomial_1.44.0 rpart_4.1.23
[137] glue_1.7.0 cachem_1.0.8
[139] Hmisc_5.1-2 generics_0.1.3
[141] Biostrings_2.70.3 classInt_0.4-10
[143] mvtnorm_1.2-4 spdep_1.3-3
[145] metagenomeSeq_1.43.0 statmod_1.5.0
[147] biomformat_1.30.0 ScaledMatrix_1.10.0
[149] carData_3.0-5 minqa_1.2.6
[151] ANCOMBC_2.4.0 glmnet_4.1-8
[153] utf8_1.2.4 seqinr_4.2-36
[155] gtools_3.9.5 readxl_1.4.3
[157] datawizard_0.10.0 ggsignif_0.6.4
[159] gridExtra_2.3 shiny_1.8.1.1
[161] GenomeInfoDbData_1.2.11 energy_1.7-11
[163] rhdf5filters_1.14.1 parameters_0.21.6
[165] RCurl_1.98-1.14 memoise_2.0.1
[167] rmarkdown_2.26 scales_1.3.0
[169] gld_2.6.6 svglite_2.1.3
[171] phylobase_0.8.12 rstudioapi_0.16.0
[173] cluster_2.1.6 hms_1.1.3
[175] ggcorrplot_0.1.4.1 munsell_0.5.1
[177] colorspace_2.1-0 rlang_1.1.4
[179] s2_1.1.6 DelayedMatrixStats_1.24.0
[181] sparseMatrixStats_1.14.0 circlize_0.4.16
[183] scuttle_1.12.0 mgcv_1.9-1
[185] xfun_0.43 multcompView_0.1-10
[187] coda_0.19-4.1 e1071_1.7-14
[189] TH.data_1.1-2 plotROC_2.3.1
[191] iterators_1.0.14 emmeans_1.10.1
[193] abind_1.4-5 tibble_3.2.1
[195] treeio_1.26.0 gmp_0.7-4
[197] Rhdf5lib_1.24.2 DECIPHER_2.30.0
[199] bitops_1.0-7 Rdpack_2.6
[201] ps_1.7.6 promises_1.3.0
[203] RSQLite_2.3.6 sandwich_3.1-0
[205] DelayedArray_0.28.0 proxy_0.4-27
[207] Rmpfr_0.9-5 compiler_4.3.3
[209] statsExpressions_1.5.4 prettyunits_1.2.0
[211] boot_1.3-30 beachmat_2.18.1
[213] BiocSingular_1.18.0 units_0.8-5
[215] MASS_7.3-60.0.1 kSamples_1.2-10
[217] progress_1.2.3 uuid_1.2-0
[219] BiocParallel_1.36.0 insight_0.19.11
[221] R6_2.5.1 rstatix_0.7.2
[223] fastmap_1.1.1 multcomp_1.4-25
[225] vipor_0.4.7 rsvd_1.0.5
[227] ggstatsplot_0.12.3 nnet_7.3-19
[229] gtable_0.3.5 KernSmooth_2.23-22
[231] latticeExtra_0.6-30 deldir_2.0-4
[233] htmltools_0.5.8.1 RcppParallel_5.1.7
[235] bit64_4.0.5 lifecycle_1.0.4
[237] processx_3.8.4 spData_2.3.0
[239] nloptr_2.0.3 callr_3.7.6
[241] vctrs_0.6.5 zeallot_0.1.0
[243] ggfun_0.1.4 sp_2.1-4
[245] pillar_1.9.0 gplots_3.1.3.1
[247] prismatic_1.1.2 ShortRead_1.60.0
[249] locfit_1.5-9.9 jsonlite_1.8.8
[251] expm_0.999-9 GetoptLong_1.0.5